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Abstract 

The network flow optimization approach is ofl'ered for restoration 
of gray-scale and color images corrupted by noise. The Ising models 
are used as a statistical background of the proposed method. The new 
multiresolution network flow minimum cut algorithm, which is espe- 
cially efficient in identification of the maximum a posteriori estimates 
of corrupted images, is presented. The algorithm is able to compute 
the MAP estimates of large size images and can be used in a concur- 
rent mode. We also consider the problem of integer minimization of 
two functions C/i(x) = XYlilVi ~ \ + Sjj ~ f^2(x) = 

Si ^iiVi - + Yji,j Pijixi - Xjf' with parameters A, Xi,fii,j > and 
vectors x = (xi,... ,Xn),y = (yi,.-. ,yn) G {0,... ,L - l}". Those 
functions constitute the energy ones for the Ising model of color and 
grayscale images. In the case L = 2 they coincide, determining the en- 
ergy function of the Ising model of binary images, and their minimiza- 
tion becomes equivalent to the network fiow minimum cut problem. 
The efficient integer minimization of Ui{x), C/2(x) by the network flow 
algorithms is described. 
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1 Introduction 



We present a new multiresolution algorithm for finding the minimum network 
flow cut and methods of efficient integer minimization of the function Ui{x.) = 

^Y.i \yi-Xi\ + J2i,j andf/2(x) = Y.i K{yi-Xi)'^ + Y.i,j P^,ji^i-^j)'^■ 

^he problem is posed in terms of Bayesian approach to image restoration 
to have unified canvas of presentation of the results, and since the results 
developed were tested and turned out efficient for processing of corrupted 
images. Though the minimization of the functions t/i(x) and Ui{x.) concerns 
not only the MAP estimation. It can be understood, for instance, as h- and 
/2-regression and so on. 

The restoration of degraded images is a branch of image processing that 
is now extensively studied for its evident practical importance as well as 
theoretical interest. There are many approaches to solution of the problem. 

We consider here the Maximum a posteriori (MAP) Bayesian estimation 
0, ^. As usual, in this statistical method one supposes an unknown image 
to be corrupted by a random noise, and that a corrupted version of the image 
is only observable. The unknown original image is presented as a random 
realization of some Markov random field (MRF), which is independent of 
the random noise. The distributions of the MRF and the random noise are 
assumed to be known. The MAP estimate is determined as the Maximum 
likelihood evaluation of the conditional distribution of the original image 
given the observed one. 

Because of use of an appropriate prior information, the Bayesian MAP 
estimators are among those, which give the best quality of restored images. 
The theory of evaluation of the MAP estimates is also significantly developed 
B S 0' il • '^^^ of disadvantages of almost all known algorithms is their 
NP-hardness (their execution time is exponent in a degree of image size). 

The first efficient method for evaluation of the Ising MAP estimate of 
binary images, which requires polynomial execution time in number of image 
pixels, was proposed by Griege, Porteous and Seheult The authors used 
the network flow optimization approach to reduce problem to identification 
of the maximum network flow. The finding of the exact Ising MAP estimate 
of binary image required 0{n^) operations, where n was an image size. But 
known maximum network flow algorithms still did not allow computation of 
the real binary images having size n = 256 x 256 or more. Therefore, in |Q an 
heuristic modification of the algorithm has been described. That algorithm 
was implemented not for entire images but for pieces of images with fixed 
values at their frontiers. 

The similar idea was used for the proposed multiresolution network flow 
minimum cut algorithm (actually, to find the MAP of a corrupted image one 
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need only the minimum network flow cut). Some additional theoretical rea- 
sons allowed describing a new algorithm, which identifies the exact minimum 
cut using special modifications of subnetworks of the partitioned original net- 
work. This algorithm can be exploited for an arbitrary network. At worst it 
takes 0{n^) operations for computation of the minimum network fiow cut, 
but for networks that correspond to real corrupted images it usually requires 
only 0{n) operations, which would be implemented in a concurrent mode. 

The description, theoretical ground and the use of the multiresolution 
network fiow minimum cut algorithm (MNFC) is presented in Section ^ 

The brief entry to two Ising models of color and grayscale images is done 
in Section |^. Those models are widely used for reconstruction of corrupted 
images, they coincide when images are binary ^, 0- In the same Section 
the problem of finding the exact MAP estimates is formulated in terms of 
the integer programming. 

In Section ^ methods of the integer minimization of functions (which are 
the energy functional of the considered Ising models) t/i(x) = A \yi — 

+ - Xj\ and f/2(x) = J2i ^iiVi - ^if + Y.i,j - XjY with 

parameters X, Pij > and vectors x = (xi, . . . ,x„),y = {yi, ■ ■ ■ ,yn) £ 
{0, ... ,L — 1}" by the network fiow algorithms (including the MNFC) is 
described. 

Applications of developed methods is given in Section |^. 

2 Multiresolution network flow minimum cut 
algorithm 

2.1 Preliminaries and Notations 

We now describe the multiresolution network flow minimum cut algorithm. 
For some types of networks this highly parallel method turns out more speedy 
and more efficient in comparison with known maximum network fiow algo- 
rithms. It was successfully used for identification of minimum cuts of large 
networks (for instance, for determination of the MAP of binary images [§) 
while classical methods were not able to solve the problem. 

The essence of the MNFC is the following. A network is partitioned 
into several subnetworks of appropriate sizes. Every subnetwork is modified 
in a special way and evaluated two times. The first time all boundary arcs 
going from the outside into the subnetwork are considered as going from 
the source (to the same nodes) and all boundary arcs going from within the 
subnetwork outside are ruled out the consideration. The minimum cut of 
the modified subnetwork is determined by known algorithm. It is proven 
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that all nodes, which are connected with the sink by a directed path, will be 

connected with the sink in the solution of the original network by the same 

path. The second time all boundary arcs going from the outside into the 

subnetwork are excluded from the consideration and all boundary arcs going 

from within the subnetwork outside are supposed to be connected with the 

sink. The minimum cut of the modified subnetwork is identified once more. 

This time all directed paths going from the source are determined. The paths 

identified are excluded from the further consideration. The reduced network 

01 is divided into several subnetworks again. The procedure of identification 

of arcs belonging to but not belonging to the minimal cut is repeated 

now with one difference - we take in account arcs connected with excluded 

nodes. The arcs not belonging to the minimal cut of 0i are found and are 

removed anew. We go to higher levels until we obtain the network &k that 

can be solved by a usual maximum network flow algorithm. 

In more details. Let the network (3 consist of + 2 numbered nodes 

S = {0,1,2 ... ,n,n + 1}, where s = is the source, t = n + 1 is the 

sink and S = {1,2... ,n} are usual nodes. The set of directed arcs is 

A = : S}. Capacities of arcs are denoted by dij > 0. Suppose 

the network satisfies the condition: 

For every usual node i & S there is either an arc (s, i) & A connecting 

(G) the source s with i or an arc {i,t) G A connecting i with the sink t 

but not both of these arcs. 

Remark 1. The condition (G) does not restrict the use of the MNFC. It has 

been done to simplify writing down notations and proofs. Any network 25 

can be easily modified to satisfy (G) for 0{n) operations (one look through 

usual network nodes S). The modified network will have the same minimum 

cut. 

Let two sets W,B d S such that W[^B = 0andW\jB = S he a 
partition of S. The vector x : S ^ {O,!}''^! with coordinates = 1, if 
i & W, and Xj = otherwise is the indicator vector of the set W. The set of 
all such indicator vectors x is denoted by X. The capacity of the s — t cut 
C(x) is defined as the sum of capacities of all forward arcs going from the 
set W U {s} to the set B U {t} ||, i.e. 

i€WU{s} jGBU{t} 

Let y be the vector with coordinates i/i = 1, if (s, i) G A, or i/i = 0, if 
{i,t) G A (remind that the condition (G) is vahd). Set Aj = dg^i, if {s,i) G A, 
or Aj = di^t, if iht) G A. The capacities of usual arcs {i,j), i,jES will be 



4 



denoted by /3ij = dij. Then in new notations 



C(x) = ^ A, + 5^ A. + 5^ A 

j e w, ieB, {i,j)esxs 
yi = o Vi = 1 



Now introduce the function 

f/(x) = ^ Ai(l - 2yi)xi + ^ - Xj)xi. 

i&S (ij)G5x5 

It is easy to see that C(x) = f/(x) + Xlj'Li -^j?/*- Since the term Yl^=i ^iVi 
does not depend on x, the functions C (x) and U (x) has minima in the same 
points X* = argminxexC(x) = argminxgxf^(x). Therefore, the solutions 
X* = argminxgxt^(x) entirely identify minimum network flow cuts p!0[ . 
For an arbitrary subset of usual nodes E G S we define two functions 

UEix) = '^\i{l -2yi)xi + ^ Pij{xi-Xj)xi (1) 

i&E {i,j)eiExE)U(ExE'')U{E''xE) 

and 

Ve(x) = ^ Ai(l - 2yi)xi + ^ Pij{xi - Xj)xi, 

ieE<= (i,j)e{E<^xE'^) 

the sum of which is equal to U (x) = f/^(x) + VE(x). Also we define restriction 
xg of X onto the set E such that x = (x^;, X£;c). The function Ve;(x) depends 
only on x^; (i.e. V£;(x) = Vj;(x£;)), therefore the following simple proposition 
is valid. 

Proposition 1. //x^ minimizes the function (t){^E) = U{xe,^e<') for some 
fixed X£;c, then for any set D C E the restriction x.'^ of x.'^ minimizes the 
function ipixf)) = t/£)(x£), xL ^, X£;c) and 



vice versa. 



For vectors Xe,Ze we will write: 
Xe < Ze, if Xi < Zi, i e E; 

xe < Ze, if Xi < Zi, i E E and there is at least one node j E E such that 



Xj < Zj] 



xe ^ze, if there are nodes i,jEE such that Xi < Zi and Xj > zj. 
Our method is based on monotone (in some sence) dependence of restrictions 
x^ of solutions x* = argminxexf^(x) on values x^c. 
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2.2 Properties of Minimum Cuts 

Let us consider the property of the monotony of x|j in details. For fixed 
and Zgc set x'^ = argminxgf/(x£;, X£;c) and = argminxgf/(x£;, Z£;c). It is 
clear that every set {x'^}o , {z'^}o can consist of more than one element. 

In general, in the case X£;c<Z£;c the inequality x^ < z'^ is not satisfied, but 
without fail there exist solutions x.'^ and z'^ such that x^ < z^. 

Theorem 2. Suppose that yiE'^ ^z^jc are fixed and let x.'^ ^ z^ (^z.e. solutions 
Xg anc? z'^ are not comparable) . Then for nonempty set D = {i a E : x[ = 
1, z[ = 0} the vector (z'^,x'^\^^) = (0£),x'^^j^) minimizes the function U 

under condition Xec, and the vector (x^,z^^^) = (l£),z^^^) minimizes the 

function U under condition Ze^, that is (z^,x^^^) G {x^}o ^, (x'^,z'^^^) G 

Proof. Suppose x^ G {x'£;}o and z'^ G {z'e}° are two incomparable solu- 
tions Xg ^ z'^ for frontier conditions xec<ze<:. It follows from Proposition |l|, 
equalities x'^ = l/j and z^ = O^, and the inequality (x'^^^, X£;c) < (z'^^^, z^;^ 
) that 

f/z)(z'^, Z^\^, Zijc) < [/^(X^, Z^\^, Zijc) < f/D(x'^, y^E\D^ ^E-) ■ (2) 

And by analogy, 

[/^(X^, X^\^^, X£;c) < Uoi^D^ ^E\D^ ^E-) < Uoiz'^, Z^\^^, Zec) (3) 

Gathering estimates in one chain we can see that all inequalities in (^,^) 
are actually equalities. This fact and Proposition |T] finish the proof. □ 

The following Corollary, which characterizes some structural properties of 
the set of solutions {x*}, is deduced from Theorem |[ 

Corollary 3. (i) If fixed frontier vectors satisfy the condition x^jc <zec, 
then for any solution x'^ = argminxj^U (xe , Xe^) there exists a solution z'^ = 
argminxj^U{xE,ZEc) such that x'^ < z'^. 

(a) If fixed frontier functions satisfy the conditionxEc <zec, then for any 
solution z'^ = argmin^j^U{xE,ZE<:) there exists a solution 
x'^ = argmin^j^U{xE,XEc) such that x'^ < z'^. 

(Hi) For any frontier condition xec the set {x'^}o has the minimal (the 
maximal) element x^^ (^e)- 

(iv) If Xec <zec, then x^ < z^ and x!^ < z'^. 
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Sentences (z) and {ii) follow immediately from Theorem]^. Sentence {Hi) 
is deduced from Theorem |^ for :kec=zec. And (iv) follows from and 
definitions of minimal and the maximal elements. 

To identify a local solution = argminxg?7(x£;, X£;c) the network 0^ 
with nodes Se = {E, {s}, {t}} and arcs i,jGS having capacities 

dij = Pij, i,jeE; ds,i = XiVi + Pj,u = Ai(l - l/Q + A,j 

1 ^ 

(4) 

is used, since the following proposition can be easily proved. 

Proposition 4. The local solution x^ sets the minimum cut of the network 
0^ and vice versa. 

The subsets E & S are chosen so that the maximum flow in the network (5'^ 
can be computed by usual maximum network flow algorithms. 

2.3 The MNFC Algorithm 

The main idea of the MNFC is to estimate at least parts of restrictions Xp of 
solutions X* = argminx?7(x) for a suitable partition f]Ei = S, Ei f]Ej = 0. 
For this purpose the monotone dependence of local solutions x^. on the 
frontier values of x*^c (in the sense of Corollary |^) is exploited. The parts 

i 

of x|j. estimated by the special local solutions x^. are ruled out the further 
consideration. It significantly reduces computational expenses. 

More precisely. Let sets Ei{l),... ,Ek^{l) partition the set 5* so that 

f]';l-^Ei{l) = and Ufii^i(l) = S. It follows from Proposition |l] that 
= argminx^.^^jf/(x£;^(i),x;^^,(^)) = x;^.^^), and the solution x^.^^^ min- 
imizes the function UE-{i)i^E-{i),^*E'=(i)) ^'^^ frontier condition x^c^-,. 
Corollary 1^ guarantees existence of two solutions (that can be found by known 
maximum network flow algorithms) 

xo,E^(i) = argminx^.(,)?7(xi5.(i),0sc(i)) 

and 

xi,i?.(i) = argminx^.(,)f/(xi5.(i), 1e9{i)) 
satisfying the inequality 

Xo,£;-(i) < x'e-{i) < ^i,E^{i) (5) 
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Denote the sets on nodes -Bi(l) = {k E -Ej(l) | xi^E-{i),k = 0} and PVi(l) = 
{k G Ei{l) I Xo,£;.(i),fe = 1}. The equahties Ob^{i) = ^i,b-{i) = x^.(^) and 
liyj(i) = Xo,Wj(i) = x^.(i) is inferred from (1). 

If the set R{1) = U£i(W^i(l) U ^i(l)) is not empty, then we identified 
the part of the solution x^^-^-j. There is a sense to continue the MNFC . 
Assign 5(2) = S\R{1) and consider now the reduced problem of identifica- 
tion x*(2) = argminxs(2)f/(x5(2),x^(^)) = argminx<,,2)t/5(2)(x5(2), x^^^^) for the 

frontier condition x^^^^. Partition S{2) = [jj^ ^^(2), ^^(2) fj ^j(2) = and 
estimate 

^'e^{2) = ^%i2) = argminxg.(2)?7E.(2)(x£;.(2),X5(2)\£;.(2),X;^(l)) 

by the solutions 

xi,£;.(2) = argminx^.,2jf/£;.(2)(xE.(2), l5(2)\i<;-(2),x;^(i)) 

and 

Xo,i?.(2) = argminx^^(2)f^Si(2)(xi?.(2)05(2)\i?i(2),x^(i)) 

satisfying the inequality xo,£;^(2) < ^'e-{2) — ^i,Ej(2)- Then, consider sets 
Bi{2) = {ke Ei{2) I xi,£;.(2),fc = 0} andVi(2) = {k e Ei{2) \ xo,E-(2),k = 1}, 
which identify x^.^2) = 0-6^(2) and x^.(.2) = 1vf^(2) correspondingly. If the set 

i?(2) = [Ji=i{Wi{2)) |Ji?j(2)) is nonempty, the algorithm is employed at the 
third level. The MNFC is iterated at higher levels until the problem will be 
completely solved or until R{1) ^ 0. At the uppermost level an appropriate 
known maximum network flow algorithm (MNF) is applied. 

Let us describe the MNFC in brief. 
Step 1 (Initialization). Assign: the level number / = 1, the set not estimated 
nodes S{1) = S, the set of estimated nodes R = 0, and the set of arcs 
A{1) = A. 

Step 2 (Partition of the network). If the maximum flow in the network 
{S{1), {s}, {t}, A{1)} can be efficiently identified by usual MNF, we do not 
need its partition. Assign the number of partition elements ki = 1, and 
Ei{l) = S{1). Identify the maximum fiow in {S (l) , {s} , {t} , A{1)} , STOP. 
Otherwise, partition the set S{1) by not intersected sets i?i(/), _E2 (/),•• • , 

Ekiil), ntiHl) = 0, UtiES) = Sil), so that every network {E^), 
{s},{t},Ai{l)} with arcs Ai{l) = Ae.{1) = {(/i,z/) : /i E Ei{l) oi u e Ei{l)} 
can be evaluated by an appropriate maximum network flow algorithm (for 
instance, by the labeling algorithm or other ones). 
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Step 3 (Computation of local estimates xo,Ej(0) Compute the vec- 

tors 

Xo,£.(/) = argminx^.(,,f/(xB.(;),05(,)\i=;.(i),x^) 

< Xi,E.(/) = argminx^,„)t/(xi5.(;), lsii)\E^(i),^*R) 

(see Corollary Q and Proposition ^) by an appropriate MNF. These two vec- 
tors correspond to minimum cut of the networks {Ei{l) , {s} , {t} , Ai{l)} for 
frontier conditions {Os(i)\Ei{i),^*R) and (l5(i)\E.(i), x;^) respectively. 
Step 4 (Identification of the set R{1) of estimated nodes). Find the sets of 
nodes Bi{l) = {i e Ei{l) \ xi^E-{i),i = 0} and Wi{l) = {i e E^{1) \ xo,E^{i),i = 1} 

keeping their values in x*. Assign the set R{1) = IJi==i(^*(0 U-^«(0) 
Step 5 (Check whether the multiresolution approach can be continued). If 
R{1) = 0, interrupt execution of the MNFC and try to identify x^^.^^ = 
argminx£,(;jf/(x5(;), x|j) by an appropriate MNF, STOP. 
Step 6 (Jump to the higher level). Assign R := R[jR{l) and S{1 + 1) = 
S{l)\R. If S{1 + 1) = 0, the proplem is solved - STOP. Otherwise, go to the 
higher level, i.e. assign / := / -|- 1, A{1) = {(yU, z/) : G S{1) or z/ G S{1)} 
and go to Step 2. 

The Step 3 of the MMC can be efficiently executed in a concurrent mode. 

2.4 The number of operations and execution time 

The number of operations required to execute the MNFC essentially depends 
on properties of the initial network 0. In the case -R(l) = the MNFC is 
reduced to the usual maximum network flow (MNF) algorithm. But there are 
a lot of networks admitting efficient application of the MNFC. Applications 
of the algorithm to the decision of real problems are described in Section |^. 
In those applications identification of network minimum cuts have required 
only 0{n) operations. 

In the following Proposition we formulate the simple necessary and suf- 
ficient condition in order to the MNFC does not degenerate to the usual 
maximum network flow algorithm. 

Proposition 5. The set Wi{l) ^ 0, respectively, Bi{l) ^ if and only if 
there exists such a set D C Ei{l) for which the inequality 

Ai(l — 2yi) + I3i^j < 0, respectively, 

i (i,j)&(DxD'^) 

5^A,(l-2|/,)+ /^M^O (6) 

i (i,j)e{D''xD) 
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is valid. 



We will write down estimates of operations number A^^op and time of 
parallel execution Tpar in general form and consider several special cases. 
Let ni{l) be number of nodes in a subnetwork, mi{l) be number of arcs in 
a modified subnetwork, and ai{l) be number of boundary arcs such that 
ai{l) = /i e Ei{l),u e Ef{l) or 1/ G Ei{l),fi G E-{1). It is sup- 

posed that the network satisfies the condition (G) (recall that modification 
of the network to satisfy this condition needs 0{n) operations). The amount 
of operations A^^op has the same order as number of operations required to 
construct all local estimates :kq^e-{i),^i,e-{i)- But the computation of each 
local estimate requires 0(aj(/)) + 0(n^(/)) operations, where 0(aj(/)) oper- 
ations are needed to modify subnetwork capacities (see, (^) and 0{n^{l)) 
operations are required directly to identify the local estimates (the second 
term 0{ni{l)^) depends on the used maximum network flow algorithm. It 
can be replaced by 0{ni{l)mf{l)) etc.). If the number of levels required 
to compute the minimum cut is equal to L, then the number of opera- 
tions is Nop = O [J2f=i Yl'iLi{(^i{^) + ^K^))) parallel execution time is 



First we consider the worst case. Suppose that all usual nodes S are com- 
pletely connected, i.e. every two usual nodes are connected by two directed 
arcs. And suppose that a prior information about a strategy of partitions 
is not known. Then, the pyramidal approach allows to avoid excessive com- 
putational expenses. At level 1 < / < logg n = L partition the original 
network into 2^^^\ k{l) = logg^i — / subnetworks containing the same number 
of nodes. Therefore, at worst, the number of operations is A"op = O(n^) and 
the parallel execution time is Tpar = 0{n^). The amount of operations will 
be O(n^) if at each level the small number of nodes is identified only, and 
the others are determined at the uppermost level. Even in this case com- 
putational expenses of the MNFC are not excessive in comparison with the 
traditional MNF algorithms. 

The situation changes if there are large subnetworks that satisfy the con- 
dition (^. The MNFC becomes very efficient. The number of operation 
can amount to 0{n). Let us, for instance, consider networks that arise from 
the Ising MAP estimation. Suppose that usual nodes 5* form the finite 2D 
lattice 5* = {1, ■ ■ ■ , A^} x {1, ■ ■ ■ , A^}, and that every usual node has 2D in- 
dex i = {ii,i2)- The nearest-neighbor nodes are connected by two differently 
directed arcs and with the same capacity Ptj = (3. Beside that 
suppose that each node are connected either with the source s or with the 
sink t by directed arcs having capacity A. Obviously, in this case the condi- 
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tion (^) will be satisfied for subnetworks with nodes having constant values in 
y and such that the set of these nodes possesses small enough isoperimetric 
ratio. Networks that correspond to corrupted binary images usually contain 
large (or even, so called, gigantic) submetworks satisfying mentioned above 
property. Those subnetworks are often estimated by the MNFC and ruled 
out the further consideration at the first level of execution of the algorithm. 
So that evaluation of the MAP estimate takes fractions of a seconds while 
the classical maximum network flow algorithms can not identify the solution 
at all. To justify this arguments we formulate the following sentence. 

Proposition 6. (i) Let By = {i e S \ yi = 0} or Wy = {i e S \ yi = 1} be 
connected component of usual nodes of the constant value in y. If P < A/2 
and By (resp., Wy) has a contour, By (resp.,Wy) maintains the same value 
in the minimum cut x* . 

(ii) If for some M > the capacities satisfy the inequality (3 < 2(m+i) ^ ^^^^ 
every connected component By (resp., Wy) of constant values iny, consisting 
of more than M nodes, maintains the same value in the minimum cut x* . 

In the case (3 > \/2 the MNFC also easily finds the MAP estimates of real 
binary images (that is equivalent to identification the minimum cut of the 
large network). The results of the practical restoration of corrupted images 
by the MNFC are placed in Section ^. 

3 Two Ising models of color and grayscale im- 
ages 

We describe two simplest Ising models of color and grayscale images with 
energy functions ?7i(x) = XY^i " + Sjj Ajki - Xj\ and t/2(x) = 
Yli^ihji - ^if + Yi,jl3i,j{xi - XjY with parameters X,Xi,l3ij > 0. Those 
models are well investigated H, 0, H and quite often applied for practical 
image restoration. The main our goal is the developing efficient algorithms 
of their integer minimization. 

Let now 5 = {1, . . . , ra} x {1, . . . , n} denote the finite square lattice of 
real numbers; then i = (^l, ^2) ^ S denotes points or pixels of our images and 
later on nodes of appropriate networks. The gray scale image is a matrix 
^gray with nouuegative integer coordinates Xi G {0, . . . , L — 1}, L > 2. The 
term "binary image" is used if coordinates of the matrix take two values 
Xi e {0, 1}. The color images Xcoi are specified by three matrices (x^, x^, x;,). 

First, we remind in brief of the Ising model of corrupted binary images 
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B B 0' HI- "^^^^ Bayesian approach specifies an a priori distribution 
p(x) = c • exp { - ^ Pijixi - XjY} 

(with the normahzing constant c and jSij > 0) over all binary images x = 
Xbj„. The conditional probability of the original binary image Xfej„ given the 
corrupted version y = is represented in the form 

p(x|y) = d{y) • exp { - ^ Xi{yi - Xif - ^ f3i,j{xi - Xjf}, (7) 

ieS i,j&S 

where the function (i(y) is independent of x. The MAP estimate, which is 
defined as a mode of the conditional distribution x* = argmaxxp(x|y), is 
equal to 

X* = argminx< ^ Xi{yi - Xif + ^ (3i^j{xi - Xjf \ (8) 

Remark 2. Often, to recover real binary images it is supposed Aj = A and 
(3ij > only for i,j belonging to local neighborhoods (under the last condi- 
tion the a prior distribution p(x) becomes the Markov Random Field). We 
knowingly consider more complicated model, since even in this case there is 
an efficient solution of the problem. 

For binary images the functions ?7i(x) and f/2(x) coincide and the sum in 
(0,|) is equal to both of them. But in the case of grayscale images f/i(x) 7^ 
t/2(x). Therefore, for x = yigray, y = ygray "we consider two different Bayesian 
models with the a posteriori probabilities 

Pi(x|y) = di(y)-exp{-[/i(x)} (9) 

and 

P2(x|y) = rf2(y)-exp{-f/2(x)}. (10) 

Both of these models are interesting from the theoretical and applied points of 
view. The first one is less investigated theoretically but it gives better quality 
of restored images. Identification of the MAP for the first model is equivalent 
to finding x* = argminx?7i(x), and for the second - Xg = argminxt/2(x), 
where argmin is taken over all grayscale images. 

For color images Xco« = (xr,Xg,X5) we consider the Ising models (|9|, p!0|) 
for each component separately (with x = x^, x = x^, x = x;, and respectively 
y = yr,y = y^^y = yt)- Thus, the problem is the same: to find x]; and Xg 
for every component of a color image. 
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4 Integer minimization of Ui{x.) and ?72(x) by 
network flow algorithms 

The network flow methods can be successfully applied for the efficient inte- 
ger minimization of Ui{x.) and U2{'x.). In particular, the MNFC, which has 
been described in Section turns out very efficient for finding of integer so- 
lutions X* = argminxgr„y?7i(x) and X2 = argmiux^^^y ?72 (x) in grayscale (and 
color) image restoration problems. The same idea is used in both cases: it 
is to represent of x with integer coordinates Xi G {0, ... ,L — 1}, L > 2 
by L — 1 matrices (for general minimization problem - by vectors) having 
0, 1-coordinates. 

4.1 Efficient minimization of ?7i(x) 

Let for integer < / < L — 1 the indicator functions l(j;->^) be equal 
to 1, if Xi > I, and be equal to otherwise. Then any grayscale image 
X is represented as the sum x = J2l'=i ^(0 of binary images-layers x(/) 
with coordinates Xi{l) = 1(2;. >/). Those binary images-layers specify the 
monotone decreasing sequence x(l) > x(2) > . . . > x(L — 1), and vice versa, 
any sequence of binary images x(l) > x(2) > . . . > x(L — 1) determines the 
grayscale image x 

= Eti We will use the following obvious 

Proposition 7. For any integer a,b & {0, . . . ,L — 1} the equality 

L-l 

k - ^1 = I l(a > ~ ^{b > I) I 
1=1 

is valid. Therefore for any grayscale images x, y the function f/i(x) can be 
written in the form 

L-l 

f/i(x) = J]nKx(/)), (11) 
1=1 

where functions 

ui{x{l)) = \^\yi{l) ~ Xi{l)\+ pij\xi{l) - Xj{l)\. 
ies («,i)eS 

Let 

x{l) = argmin3;^.^u;(xbin), / G {1, . . . , L - 1} (12) 

denote binary solutions that minimize the functions 'Ui(a;bin)- Since y(l) > 
y(2) > . . . > y{L — 1) the following sentence is valid. 
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Proposition 8. There is a monotone decreasing sequence xm(1) > ^m(2) > 
. . . > xm{L — 1) of solutions of (fZ^j. 

The proof of Proposition ^ is similar to the proof of Theorem ^ in Sec- 
tion ^ To show existence of solutions x]; of the form x]; = x.1 = J2i=i ^m{1) 
it is enough to use Propositions ^ and |^. Indeed, for any x 

L-l L-l 

1=1 1=1 

Each binary solution xm{1) can be identified by the maximum network flow 
algorithms or by the MNFC for polynomial number of operations. At worst, 
it takes 0{Ln^) operations but quite often the use of MNFC allows reducing 
operations up to 0{Ln) operations that are held in a concurrent mode (see 
Section 

4.2 Efficient minimization of t/2(x) 

The main idea of the efficient integer minimization of the function f/2(x) = 
Aj(?/j — XiY + J2i j i^iji^i ~ ^jY li^s replacement of the variable x by 
the sum of unordered Boolean variables x(/), and then considering another 
polynomial g(x(l), . . . ,x(L - 1)) > f/2(Eti ^(O) of 

many Boolean vari- 
ables with points of minimum q(l), . . . , q(L — 1) such that Xg = Ylt=i ^(0 
minimizes t/2(x). The new polynomial Q is chosen so that it admits the 
efficient minimization by the network flow optimization algorithms. 

In more details. Let us represent x as the sum of arbitrary Boolean 
variables x = Ylt=i ^(0- ^^"^ variables the polynomial t/2(x) will look as 

[/2(X) = E A. (y^ - E^^(O) + fE(^^(0 - 

i \ 1=1 J i,j \l=l 

Assign for simplicity bij = and use the equality x'f{l) = Xiil) for Boolean 
variables Xi{l) to write t/2(x) in the form 

f/2(x) = E >^^y^ + . . . , X(L - 1)), 

where the polynomial of many Boolean variables 
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P(x) = P(x(l),... ,x(L-l)) 



= E 



ies 



Xi - 2Xiyi - (L - 2) + 



L-l 



1=1 



i£S L ieS J l</<m<L-l 

L-l . ^2 

+ XI ^mXI 

ije5 1=1 ^ 



^ a;i(/)a;i(m) 



(13) 



ij65 l<i<m<L-l 



{xi{m) - Xj{l)y + {xi{l) - Xj{m)f 



Intoduce another polynomial of many Boolean variables 

g(x(i), ... ,x(L-i)) 

Xi - 2Xiyi - (L - 2) Y^^hj + h 



L-l 



Z=l 



je5 L jG5 J l<l<m<L-l 

L-l / X 2 

i,jes 1=1 ^ 



(14) 



ije5 l<;<m<L-l 



{xi{m) - Xj{l)Y + (xj(Z) - Xj{m)f 



such that (5(x(l), . . . ,x(L — 1)) > P(x) and which differs from P(x) by the 
term X]i<i<m<L-i ™ second summand. 

Denote by 

(q*(l),q*(2),... ,q*(L- 1)) = argmin^(i),^(2)„„,^(L_i)(5(x(l),... ,x(L- 1)) 



any collection of Boolean vectors that minimize (5(x(l), . . . ,x(L — 1)). We 
will prove that without fail q*(l) > q*(2) > ... > q*(L - 1). This feature 
will allow expressing solutions of the initial problem as X2 = Y^]^=i q*(0- 
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Theorem 9. Any collection (q*(l), q*(2), . . . , q*(L — 1)) that minimizes Q 
forms the decreasing sequence. It identifies the solution of the original prob- 
lem by the formula X2 = J2f=i '^''^^ vise versa, each solutions spec- 
ifies the solution q*(0 = x*(/). 

Proof. Let us represent Q as 



g(x(i),... ,x(L-i)) = 



^ (l - Xi{l))xi{m) 



ieS L jeS J i<i<m<L-i 

where x = X^^i^ x(/) . Then suppose that there exists an unordered collection 
(q*(l), q*(2), . . . , q*{L — 1)) minimizing Q. For this collection 

Y {i-q*m:{m)>o 

l<l<m<L-l 

and, therefore, for x = J2f=i the inequality Q (q*(l), ... , q*(L — 1)) > 

P(x) holds (remind that all Xi, Pij > 0). But, evidently, for ordered Boolean 
matrices-layers x(l) > . . . > x(L — 1) with coordinates Xi{l) = > we 
have 

Q(x(l),... ,x(L-l)) = P(x). 

Hence, unordered sequence (q*(l), q*(2), . . . ,q*(L — 1)) can not minimize 
Q. Since for any ordered collection (x(l) > . . . > x(L — 1)) the equality 
(5(x(l), . . . ,x(L — 1)) = P(x) is valid, actually, X2 = J2f=i ^ise 
versa, each solutions Xg specifies the solution q*(/) = x*(/). □ 

It follows from Theorem ^ that P and Q have the same ordered set of 
Boolean matrices (x*(l), x*(2), ... , x*(L — 1)) that minimize these polynomi- 
als. But the problem of minimization of Q is known in discrete optimization 
|]10| , 0, 0]. It is equivalent to identification of the minimum network flow cut. 
To describe the corresponding network let us re-arrange the polynomial to 
the form 



Q(x(l) 



,x(L-l)) 



L-l 



EE 

i£S 1=1 



{21 - l)\ - 2\y, + {21 - L) Y,{K, + hi 



Xi{l) 



+ 



+ 



i,jeS 1=1 ^ ^ 



E'y E 

hjeS l<l<m<L-l 



{xi{m) - Xj{l)f + {Xi{l) - Xj{m)f 
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Then denote for brevity 



{21 - - 2Kyi + {21 - L) Y,{h,j + 6,- i) 

3&S 



and introduce the network N{s, t, V, A) with the source s, the sink t and usual 
nodes that are labeled by ii G V, where the multiindex i & S and / shows a 
layer number (there arc \V\ = {L — 1)\S\ usual nodes in the N{s.t.VjA)). 
For instance, means the node i at the third layer. The set of arcs with 
corresponding capacities is specified as 



{iljm) 



—di,i 
di,i 



if 
if 
if 



di,i < 0, 
di,i > 0, 



for / g{1,...,L — 1}. In general, it is necessary 0{Ln^) operations to find a 
solution X2, but use the MNFC or other decomposition methods can reduce 
the operation number. 



5 Applications 

It was mentioned in Introduction the MNFC turned out very efficient in 

recovering of images (both binary and grayscale) corrupted by random noise. 
In Fig. 1 the binary 512 x 512 image corrupted by rather strong Bernoullian 
noise with the parameter p = 0.3 is placed. The Ising model with the energy 
function C/i(xbin) = ?^2(xbin) was used for image restoration. 
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Fig.l Fig.2 

The result of implementation of the MNFC at the first level is drawn in 
Fig. 2. The pixels that were not estimated at the first level are depicted 
in gray color. They form small sparse sets at boundaries of squares of the 
partition. The identification of x^j^^ of the corrupted image takes fractions of 
a second while the known maximum network flow algorithms were not able 
to evaluate the estimate at all. 

Now consider the application of the proposed method of the integer mini- 
mization for restoration of corrupted grayscale images. The Ising model with 
the energy function ^7i(x) was used to construct the MAP estimates. This 
model was preferred since, flrst, the function [/i(x) admits more efficient 
minimization in comparison with f/2(x) and, second, because of better visual 
quality of the estimate. Beside the MAP estimate x^, the moving average, 
the moving median estimates and the gradient estimate of the U2{x.) were 
computed to compare estimators. 




(e) (f) 
Fig. 3 
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The original grayscale image in Fig. 3(a) was corrupted by the exponential 
noise Fig. 3(b). The results of the 3x3 moving average and 3x3 moving 
median filtration are depicted in Fig. 3(c) and Fig. 3(d) respectively. The 
continuous gradient estimate for X2 is placed in Fig. 3(e). Fig. 2(f) is the 
exact estimate of the image in Fig. 3(b). Remind that it is determined 
separately for every binary layer of a grayscale image. Therefore, the use of 
MNFC allows computation of in a highly concurrent mode. 
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